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ABSTRACT 

We present a three-dimensional, fully parallelized, efficient implementation of ionizing 
UV radiation for smoothed particle hydrodynamics (SPH) including self-gravity. Our 
method is based on the SPH/tree code VINE. We therefore call it iVINE (for Ioniza- 
tion + VINE). This approach allows detailed high-resolution studies of the effects of 
ionizing radiation from e.g. young massive stars on their turbulent parental molecular 
clouds. In this paper we describe the concept and the numerical implementation of 
the radiative transfer for a plane-parallel geometry and we discuss several test cases 
demonstrating the efficiency and accuracy of the new method. As a first application, 
we study the radiatively driven implosion of marginally stable molecular clouds at var- 
ious distances of a strong UV source and show that they are driven into gravitational 
collapse. The resulting cores are very compact and dense exactly as it is observed in 
clustered environments. Our simulations indicate that the time of triggered collapse 
depends on the distance of the core from the UV source. Clouds closer to the source 
collapse several 10 5 years earlier than more distant clouds. This effect can explain the 
observed age spread in OB associations where stars closer to the source are found to 
be younger. We discuss possible uncertainties in the observational derivation of shock 
front velocities due to early stripping of proto-stellar envelopes by ionizing radiation. 
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1 INTRODUCTION 

As hydrodynamical simulations become more and more ad- 
vanced one of the key issues is the successful implementation 
of additional physics like the effects of radiation. Prominent 
applications are for example the reionizat ion of the early 
universe (for a comparison of methods see llriev et all [20061 
and references therein). 

In our present day universe ionizing radiation still plays 
a vital role. UV-radiation from massive, young stars ionizes 
their surrounding. The hot, ionized gas then expands into 
the cold, neutral gas and thus drives shock fronts into the 
parental molecular clouds. Up to now it is not fully under- 
stood if this violent feedba c k enh ances or hinders star for- 
mation. lElmegreen fc Ladal l| 19771 ) proposed that the shock 
front builds up dense regions by sweeping up the cold gas, 
which then eventually collapse due to gravitational insta- 
bility and form stars. This is called the 'collec t and collapse 
model' (see also the review by |Elmegreen|[l998l ) . Another sit- 
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uation arises when preexisting, dense structures (e.g. molec- 
ular cloud cores) that are gravitationally marginally stable 
get compressed by the approaching front and start collaps- 
ing. TjusJs£OjruTi£rily_called radiation driven implosion (see 
e.g. ISandford et "aHll982h . 

O bservations provide widespread evidence for these sce- 
narios (|Sugitani et al . 1989, Hester et al.lll996T ). More recent 
observations indicate triggered star for mation on the edge s 
of HH-regions e.g. i n the Orion clouds (IStanke et alj 2002). 
the C arina nebula ilSmith et al.ll200ol ) M16 (iFukuda etaj] 
2002), M17 (I Jiang et al.ll2002l). 30 Dor |Walborn et al.|[2 002) 
and t he SMC (|Gouliermis et all 120071 ). iDeharveng et all 
(2005) report triggered star formation in samples of more 
distant HII regions. Besides these quite complex large scale 
regions there have been numerous observations of bright 
rimmed cometary globules. These are small isolated clouds 
with a clear head to tail structure with the dense head s 
pointing towards an ionizing source (ISugitani et al.lll99ll ). 
Their morphology enables a direct comparison to simula- 
tions. In particular, the properties of individual young stellar 
objects (YSOs) surrounding OB-associations can be deter- 
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mined precisely. YSOs are observed in the mass range from 
T Tauri (0.1 - 3M B ) up to Herbig Ae Be (2 - 8 M ) stars 
(see e.g. lLee fc Cnenll2007l . ISnider et al.ll2007h . The velocity 
of the shock front triggering the star formation is calculated 
from the age difference of the stars and their relative dis- 
tance. These estimate s are in the range of a few km/s (e.g. 
iThompson et al.ll2004 iGetman et "aill2007h . 

Numerous simulations on the topic of cloud evapora- 
tion and se q uentia l star formation have been performed. 
lYorke et al.l l| 19821 ') (and references therein) published a 
series of two-dimensional simulations on the gas dynam- 
ics of HH-regions, especially on champagne flows, where a 
stream of hot gas break s through the border of cold, con- 
fining gas. Subsequently, lElmegreen et al] (|l995l) presented 
two-dimensional, grid-based simulations showing that the 
expansion of an HH-regi on into the surrounding cloud 
can trigger star formation. iKessel-Devnet fc Burkerti (2003) 
demonstrated with a three-dimensional SPH code, that a 
marginally stable molecular cloud core can be triggered 
into collapse when exposed to strong UV radiation. With a 
more deta iled descript i on of radiation implemented into an 
SPH code lMiao et al.l |2006t ) could reproduce the observed 
features of the Eagle Nebula, including the photodissoci- 
ation regions and the temperat ure profile. Us i ng a three- 
dimensional grid-based scheme, iMellema et alj (|2006l ) sim- 
ulated the HH-region excavated by a point source of UV- 
radiation. They find remarkably similar morphologies and 
physical properties when comparing their models to ob- 
servations. Furth ermore simulations with an SPH-c o de by 
iDale et all (|2005h and a grid code bv lMac Low et all (|2007T ) 
showed that a turbulent interstellar medium surrounding an 
O-star allows the ionizing radiation to efficiently expel most 
of the nearby gas. Only the denser regions resist and con- 
tinue to collapse. 

However none of the authors so far described ionizing 
radiation as an efficien t trigger for star formation. There is 
only weak evidence bv lDale et al.l (|2007h . that the external 
irradiation of a collapsing cloud by a point source can indeed 
increase the star formation efficiency from 3% to 4% when 
compared to a control run without radiati on. For a review o f 
feedback processes we refer the reader to iMac Low! (|2007t ) . 
For completeness we would like to refer to recently published 
implementations f or ion izing radiation into an SPH code by 
IPawlik fc Schavel (|2008l ). where the photons of a source are 
followed along cones, and lAltav et al.l l|2008l ). where the ra- 
diation is followed via a Monte Carlo ray-tracing scheme 

All these studies demonstrate that there is a strong 
connection between the UV-radiative feedback from massive 
stars and the observed morphologies of the ambient molec- 
ular cloud gas. Yet, a quantitative discussion of the interac- 
tion between UV-radiation and turbulent molecular clouds 
is still missing. To advance our understanding, we introduce 
iVINE, the fully parallel implement ation of UV-radiation 
in the parallel tree -SPH-code VINE (|Wetzstein et al.ll2008l. 
I Nelson et al.1120081 ). This efficient tool permits high resolu- 
tion simulations of molecular clouds in the vicinity of strong 
UV sources such as an O-star or association. 

The paper is structured as follows. The physical model 
and its implementation are described in Section [2] followed 
by a detailed comparison of the scheme with analytical re- 
sults (Sec. [3]). We apply the new method to the radiatively 
driven implosion of a marginally stable molecular cloud core 



and compare three simulations with different initial UV 
fluxes to observations (Sec [4| . In Section [5] we summarize 
and discuss the results. 



2 NUMERICAL METHOD 

As soon as a young massive star emits UV-radiation it ion- 
izes its surrounding, creating a HH-region. Initially the ion- 
ization proceeds fast with a speed of this rarefied (or PL- 
type) front of v-r >> dhot, where ahot is the sound speed 
of the hot, ionized gas. After a sound crossing timescale the 
hot gas reacts to its increased temperature and an isother- 
mal shock front is driven into the cold surrounding medium. 
This dense (or D-type) shock travels at a much smaller speed 
vd ~ ahot. For a full t extbook analysis of this evolution see 
e.g. lOsterbrockl |l989h . 



2.1 Prescription of ionizing radiation 

To follow the evolution of the HH-region of a young mas- 
sive star in a numerical simulation we use a prescrip- 
tion for the ionizing UV-radiation simila r to the one pro- 
pose d by [K csscl-Dcynct & Burkert] (|2000f ) as presented be- 
fore (|Gritschneder et al.ll2007h . The flux J at any given po- 
sition x is given by 



J(x) = Jl y 



(1) 



where Jl y is the Lyman continuum flux of the hot star. The 
optical depth t„ is given by the integral along the line of 
sight between the source of radiation and the position x 



n v pdx, 



(2) 



where re„ is the frequency weighted absorption coefficient 



(3) 



with riH being the number density of neutral hydrogen and 
p the mass density of the gas. We assume the gas is pure 
hydrogen with a mean molecular weight of /J, = 1. As the 
frequency dependent absorption cross section a v peaks at 
the Lyman break it is a valid assumption to take an average 
cross section <r, thereby approximating the radiation to be 
monochromatic. Thus, every photon above the Lyman break 
is assumed to ionize one hydrogen atom. 
We define the ionization degree i) as 

(4) 

n 

where n c is the number density of electrons and n is the 
combined number density of protons and neutral hydrogen 
atoms. The time derivative of the ionization degree can be 
written as 

d-q _ 1 dn c 
dt n dt n 

with the ionization rate / given as 



-(I-R), 



(•») 



/ = VJ 

and the recombination rate R as 



R = n^QB = rj 2 n 2 aB- 



(6) 



(7) 
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For the recombination coefficient aa we choose 



(8) 



where ai is the recombination probability for a level i of the 
hydrogen atom. The recombination of electrons and protons 
leads to a diffuse field of Lyman continuum photons, which 
in turn can again ionize a hydrogen atom. We neglect this ef- 
fect under the assumption that every reemitted photon is in 
turn immediately absorbed in the direct surrounding. This 
assumption, called 'on the spot approximation' is valid as 
long as the hydrogen density is high enough (e.g. ISpitzerl 
1 19781 ). which is always true in the vicinity of the ioniza- 
tion front. Some fraction of the UV- photons is absorbed by 
dust, and re-emitted in the IR-regime, leading to an effective 
lower flux. We neglect this effect, since the flux incident on 
the simulation volume is determined largely by its distance 
from the radiation source, so that geometrical dilution of the 
radiation field is likely to be more important than absorption 
by dust. 

The average temperature of the gas is coupled linearly 
to the ionization degree 77 through 



0) 



where T co id is the initial temperature of the cold, unionized 
gas and That is the average temperature of the ionized gas. 

2.2 Implementation 

To treat the hydrodynamical and gravitational evolution of 
the gas we use the parallel smoothed parti cle hydrodynam- 
ics (S PH) co de called VINE dev eloped by IWetzstein et al.l 
(|2008h and iNelson et all l|2008l ). SPH is a Lagrangian 
method, which renders it extremely suitable to cover several 
orders of magnitude in density and spatial scale, for exam- 
ple during cloud core collapse and star formation. VINE is 
a powerful parallel implementation of the SPH method in 
combination with a tree code for the calculation of gravita- 
tional forces. It offers a Runge-Kutta integrator as well as 
a Leapfrog integrator. Both schemes can be used in com- 
bination with individual particle time-steps. For this work 
the Leapfrog integrator is chosen. Every time the equation 
of state is computed we calculate the ionization degree for 
all particles in the entire simulation. 

The heating by UV-radiation can be treated as decou- 
pled from the dynamic evolution since the recombination 
timescale 



ire 



J1QB 



(10) 



is much shorter than any hydrodynamical timescale. In our 
simulations (see Section 2} the crossing time even in the hot 
gas is £hot ~ 70kyr, whereas the recombination timescale is 
t TCC ~ lkyr for a number density of n = 100cm" 3 . Thus, it is 
valid to treat ionization and hydrodynamics as two separate 
processes. In other words the ionization can be assumed to 
happen instantaneously. Frequently updating the ionization 
degree together with a modified time-step criterion (see Sec- 
tion 12. 3[) ensures that the radiation is treated correctly on 
all scales. 

To include the effect of UV-radiation we assume a plane- 
parallel geometry, i.e. parallel rays. This is valid as long as 



the distance from the source of radiation is larger than the 
dimensions of the area of infall. In our simulations the radi- 
ation is impinging from the left hand side, that is from the 
negative x-direction. To couple ionization to hydrodynam- 
ics we use a ray-shooting algorithm. As the ionizing flux is 
propagated along the x-direction, we ensure the conserva- 
tion of flux by dividing the (y, z)-plane into sub-domains 
of equal size, whose extent along the ^-direction spans the 
whole simulation domain. Along each of these sub-domains 
or rays the flux is transported in a conservative manner. To 
convert the SPH-particle density p par t correctly into a gas 
density distribution within these three-dimensional rays the 
volume and thus the diameter d pa rt that each SPH-particle 
occupies is calculated via the mass of each particle m part : 
1/3 



2 ■ 



3 m 



part 



47T p pa 



(11) 



The width Ay and the height Az of the rays or bins is then 
set to the average value d pa rt of the particles closest to the 
source. To determine this value the first two particles in each 
ray at the previous time-step are taken into account. Since 
this is the region with the lowest density throughout the 
entire simulation this guarantees that the bin-size is always 
larger than the characteristic particle resolution. As soon as 
the ray approaches a density increase the local d par t becomes 
smaller than Ay and Az. For 



d r , 



cL 



1 

<2 



Ay Az 2 

we refine the ray subsequently into four sub-rays to treat the 
ionization of high density regions correctly. Each of the sub- 
rays inherits the optical depth of the main ray. Currently the 
code allows for five levels of refinement, thus increasing the 
effective bin resolution in each direction by a factor of 32. 
In principle it would be possible to de-refine the sub-rays by 
using the average optical depth of the four refined sub-rays 
for the de-refined bin. We do not include this, since it would 
lead to an unphysical shading of lower density sub-rays as 
soon as they are combined with a high density sub-ray due 
to an overestimation of the optical depth. 

To calculate the optical depth, we sort the particles 
within each bin according to their distance to the source 
and discretize into subsections of the size 



Axi 



(13) 



Thus, Axi is the projected distance of a particle i to its 
direct neighbours closer and further away from the source, 
i.e. the length along the line of sight the particle occupies. 
We then calculate the optical depth r along each ray by 
summing up the individual optical depths Tj of each particle 
i. The discrete value of Ti is given according to Eq. [5] as 



n = a uha Axi 



(14) 



The number density «h and the density p used to calculate 
the recombination rate (cf Eq. [7| are simply given by the 
SPH-density p par t. From these quantities the new ionization 
degree r\ is determined according to Eq. [5] by a Newton- 
Raphson iteration scheme. It converges with a precision of 
more than 0.1% in less than 200 iterations. When the ion- 
ization degree reaches a value of 77 = 1 x 10 -10 we terminate 
the further calculation of this bin. This implementation is 
fully parallelized in OpenMP. 
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2.3 Modification of the time-step criterion 

A detailed discussion of the different time-step criteria 
implemented in the underlying VINE code is given in 
IWetzstein et alj (|2008T ). Note that our implementation of 
ionizing radiation is designed to b e used in connection w ith 
individual particle time-steps (see IWetzstein et alJl2008L for 
details). To exactly follow the evolution of a particle during 
its ionization process it is vital to use a small enough time- 
step. To do so we decided to force every particle to a smaller 
time-step as soon as its ionization degree reaches n > 10~ 3 , 
i.e. when the particle is going to be ionized. The new time- 
step is chosen by a modified Courant-Friedrichs-Lewy (CFL) 
condition according to 



At ne 



acoid/ahot * Aic 



(15) 



where a co id and dhot are the fixed respective sound-speeds 
of the cold and the ionized gas at T co id and That- A/;cfl 
is the individual time-step the particle would get as signed 
due to the CFL-condition (see I Wetzstein et alJl200i ). This 
ensures that the hydrodynamical quantities are treated cor- 
rectly even though the particle gets a boost in tempera- 
ture. Therefore, we anticipate the subsequent acceleration 
due to the approaching ionization front by choosing already 
the much smaller time-step even though the particle is just 
ionized to 0.1%. This criterion also ensures that the ioniza- 
tion degree is followed accurately during the evolution of 
the later dense or D-type ionization front, because vu is al- 
ways smaller than the sound speed of the hot gas vo < ai wt . 
Hence, this front will always be resolved by particles which 
have a small enough time-step to track the hot gas. In the be- 
ginning the evolution of the faster R-type front (dr >> ahot) 
can be followed by using a small enough initial time-step 
since this phase is quite short (~ 5kyr). 

The choice of a small initial time-step together with the 
modified CFL-criterion ensure that the ionization degree rj 
of a particle never changes by more than ±0.1 per time-step 
in all of our simulations. Thus, the ionization front can be 
followed in both stages (R- and D-type) precisely. 



3 NUMERICAL TESTS 

In order to validate the algorithm we perform several tests. 
The first series of simulations addresses the evolution of the 
Stromgren solution and tests whether the time-dependent 
UV-flux is treated correctly on all scales. In addition, we 
demonstrate the correct implementation of the refinement 
(Section 13. ip . The second series of simulations (Section l3,2p 
is designed to demonstrate the correct interaction of ioniz- 
ing radiation and hydrodynamics. In the end the successful 
parallelization of the code is shown (Section 13. 3p . 

3.1 Ionization without hydrodynamics 

3.1.1 The Stromgren test - Ionization by a constant UV 
flux 

When hydrodynamics is not taken into account, the ho- 
mogeneous surrounding of an ionizing source will always 
converge towards an equilibrium between ionization and re- 
combination. Th e volume of the ionized Stromgren sphere 
jStromgrenll 19391 ) around an O-star is given by 



(16) 



assuming a monochromatic source with a constant UV-flux 
Jhy given in photons per second, cub and n are again the re- 
combination coefficie nt and the number density (for a text- 
book analysis see e.g. IShulll99ll ). 

In the case of plane-parallel radiation, as discussed here, 
this volume is characterized by the length x s , which can be 
penetrated by the ionizing radiation. x s is determined by 
the surface S on which the photon flux per area and time, 
F Ly , is impinging: 



Fj 



L.v 



s . . (17) 

The time evolution of the length xi(t) of this region is given 
by the differential equation 

dxi 
~dt 

with the solution 



-n — Fl y — xi(t)a^n 2 



(18) 



Xl(t) = x s (l 



-t/t r . 



(19) 



where t ICC = 1/(tiqb) is the recombination or Stromgren 
timescale. The shape of the front is given by the ionization 
equilibrium equation 



FvOvdv = n 2 ri 2 aB, 



(20) 



which can be rewritten for the plane-parallel, monochro- 
matic case in terms of the ionization degree (cf Eq. [4} as 

— = rj — — nax B . (21) 
ax 1 + 77 

This equation can be solved numerically and gives the ion- 
ization degree 77 at any given position x for the chosen num- 
ber density n and mean cross-section a. 

To test the code, we ran three simulations: 

• Case A: 125k particles placed on a Cartesian grid 

• Case B: 100k particles placed randomly 

• Case C: 250k particles placed randomly 

For cases B and C the particles are placed randomly in 
the simulation box and then are allowed to relax with peri- 
odic boundaries and the inclusion of hydrodynamics for one 
crossing timescale to dampen the numeric random noise. 
Thereafter we switch off the hydrodynamics and compute 
the ionization. For all simulations we used a mean den- 
sity n — 10cm -3 . The simulated volume is (2pc) 3 , the 
length the ionization can penetrate is set to x s = lpc. The 
recombination coefficient and the absorption cross-section 
are set to typical values of «b = 2.7 x 10~ 13 cm 3 s -1 and 
a — 3.52 x 10~ 18 cm 2 . For the above parameters, Fi,y — 
8.33 x 10 7 photons cm s and t ICC = 11.8kyr. The simu- 
lations run up to t = 5/Jrcc to allow for a quasi-equilibrium 
state to evolve. 

In Fig. [1] the time evolution of the penetration length 
xi(t) is shown. The position of the front is calculated by 
projection of the three-dimensional simulation along the y- 
and z-axis onto the x-axis. Note that the analytical solu- 
tion (cf Eq. I19p is based on the idealized assumption that 
the medium is fully ionized (n = 1.0). However, the precise 
solution of Eq. [5] in equilibrium (drj/dt — 0) is 
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t/t. 



Figure 1. Time evolution of the ionization degree n for the three 
test cases with different particle numbers and distributions red 
(125k particles), green (100k particles) and blue (250k particles) 
lines. The black solid line denotes the exact solution, the black 
dashed line the classic (Stromgren) solution. 




Figure 2. Ionization degree n (as 1 at x/x B = 0) and neutral gas 
fraction \ = 1 — n versus position for the cases A (red), B (green) 
and C (blue). The dashed line represents the exact solution. 



2 V s F Ly 

= ~o = 2 ' 



(22) 



In our simulations x B rj 2 = lpc is realized with x B = 1.05pc 
and rj — 0.976. We call this the exact solution whereas the 
solution assuming rj = 1.0 will be referred to as classic solu- 
tion. Our simulations converge very well towards the exact 
solution. Case A, where the particles are initially placed on 
a grid, slightly overestimates the final value of x$, while the 
low resolution simulation (case B) underestimates it. Never- 
theless, already with only 100k this implementation shows a 
very good agreement with the analytical curve. In the high 
resolution simulation (case C) the numerical result lies right 
on top of the predicted line. 

Fig. [2]shows the ionized fraction r\ and the neutral frac- 
tion x — 1 — V after t = 5t B at the end of the simulation. 
The numerical solution of Eq. [5T] is evaluated for the exact 
solution with a penetration length of x B = 1.05pc. As ex- 
pected from Fig. [T] case A overestimates the front position, 



whereas case B underestimates it. Again the high resolution 
run C shows the best concordance and we can conclude that 
these results fit well within the range of the code compar- 
isons done bv llliev et all l|2006h . A more direct comparison 
to this work is not possible due to the plane-parallel nature 
of the test performed here. Q 



3.1.2 Ionization by a time-varying source 

A more challenging test is the treatment of a time-varying 
source of ionization. Although this situation is not very re- 
alistic for an O-star it nevertheless provides a very good 
method to test the treatment of a rapidly changing flux by 
the code. 

To produce an ionization front that is traveling at a 
constant speed through a medium of constant density it is 
sufficient to increase the flux per area FLy linearly with time, 



-PLy(i) = UVf + n 2 CtBV[t, 



(23) 



where Vf denotes the speed of the ionizing front. The first 
term on the right hand side provides the ionization of the 
front, while the second term compensates for the loss of flux 
due to recombinations on the way towards the front. We 
assume a constant density of n — 10cm -3 and the veloc- 
ity of the front is set at vt = 1.3 x 10 5 cm s -1 . The other 
parameters are chosen as before. 

Again the three initial conditions A, B and C from sec- 
tion 13.1.11 are explored. The results are shown in Fig. [3] 
As before the simulations match the theoretical solutions 
closely. In the beginning run A agrees very well with the 
solution. This is due to the very low numerical noise in the 
Cartesian grid. However, towards the end the low resolution 
leads to a deviation from the analytic value. In case B one 
can clearly see the effect of the noisy density distribution, 
since for the recombination R any error in the density leads 
to a quadratic error in the absorption of the photons (cf Eq. 
[7]). Therefore, the position starts to oscillate around the ex- 
act position. This effect gets stronger the further the front 
penetrates, as more material has to be kept from recombin- 
ing. Case C shows a very good agreement with the exact 
solution, the resolution is high enough to keep the noise in 
the density distribution low and thus the position of the 
front is followed precisely. 



3.1.3 Testing the refinement - Ionization by a constant 
source in a two-density medium 

All tests up to now were independent of the implementation 
of refinement, since in a constant density medium each par- 
ticle occupies roughly the same diameter d par t (see Section 
I2.2|l . To verify the correct implementation of the refinement 
we set up a simulation with a two-density medium. A lower 
density gas phase with m = 10cm -3 is set up in the left half 
of the box and a higher density medium with ni = 100cm -3 



1 Note that in Fig. [2] the neutral fraction \ converges towards 
a value of 10 -2 at x = Opc in both the simulations and the 

exact solution whereas in llliev et ah | h I X is reaching much 

lower values. This is due to the fact that in our simulations the 
irradiated surface stays constant whereas when simulating a point 
source this surface and thus \ can g e * infinitesimally small. 
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Figure 3. Numerical simulation of an ionization front that moves 
with constant speed through a medium of constant density. Plot- 
ted is the position of the front in units of the box length versus 
the time in units of the crossing time for the three cases A (red), 
B (green) and C (blue). 




0.2 

x [pc] 



Figure 4. Effect of the refinement on a diagonal density step. 
Plotted are the SPH-particles in a (0.5pc) 3 volume projected 
along the z-axis. Red: ionized particles (n > 0.1), black: unionized 
particles. Left: without the inclusion of refinement. Right: with 
the inclusion of one level of refinement. 



is placed at the right half of the box. The density contrast 
is achieved via a different number of SPH-particles in the 
different regions, the particle masses are equal in the en- 
tire simulation. The required flux to ionize the simulation 
domain up to a position :r s can be calculated by linear su- 
perposition according to Eq. [33] 



2 



(24) 



where Xi = 0.5pc is the extent of the low-density region. The 
simulation is set up with 550k randomly placed particles. 
The particle noise is reduced for both regions separately as 
described in Section [3.1.11 The penetration depth is set to 
x s — lpc. As soon as the equilibrium state is reached the 
numerically calculated penetration lengths are 



t'cq, unrefined 



0.985 x s 



2-eq, refined 



0.997 x s 



(25) 



for the unrefined and the refined case. Very good agreement 
even with the unrefined code can be expected, as we always 
use the SPH-density in the calculation which is independent 
of the chosen bins for the calculation of the ionization and 
the recombination. 



Nevertheless, the refinement has an important geomet- 
ric effect, which becomes clear when assuming a density con- 
trast with a discontinuity which is not aligned vertical to the 
impinging radiation. We perform a test with a diagonal den- 
sity contrast between two regions with a number density of 
"low = 10cm -3 and JT-jjjgjj = 200cm -3 respectively. Again 
the particles are placed randomly and the noise is reduced 
(see Section l"3.1.1[l . A cubic domain of (0.5pc) 3 including 25k 
particles is shown in Fig. [4] In the unrefined case (left hand 
side) the effect of the original bin-size of « 0.05pc can be 
clearly seen as step-like features. With refinement the den- 
sity contrast of 20 leads to one level of refinement (since 
dp ar t/Ay w 1/2.7, cf Eq, [T2)l and the geometrical bias is al- 
ready negligible (Fig. [3] right hand side) . In the simulations 
in Section [4] all five levels of refinement lead to spatial reso- 
lution of the radiation in our simulations as high as 10 -3 pc, 
therefore the radiation does not produce any unphysical ge- 
ometrical effects. 



3.2 Ionization with hydrodynamics 

3.2.1 Steady propagation of an ionizing front 

This test was originally proposed by iLefloch fc Lazarefll 
1 19941 ). An area of constant density is ionized by a photon 
flux which increases linearly with time. This leads to a hy- 
drodynamical shock wave traveling at a constant speed. The 
number densities ni,n c ,no in the ionized, the compressed 
and the undisturbed medium can then be calculated from 
the corresponding sound speeds ai, a c , ao. Let iij be the speed 
of the ionization front and u s be the speed of the shock 
front. The jump condition for a D-type ionization front can 
be written as 



of, 



(26) 



since the compressed and the neutral medium have the same 
temperature and thus the same sound-speed. At the isother- 
mal shock the jump conditions are 



u s (u B — Ul) 



2 

ao, 



(27) 

- = 4. ( 28 ) 
no 

where Ui is the gas velocity just inside the shock. With the 
approximation of a thin shock it can be assumed that the 
ionization front and the shock front have the same speed 
m = u s . For a detailed derivation of the jump conditions see 
e.g. lShul (|l99ll ). Introducing the time derivative of the ioniz- 
ing flux C = &F/&t the speed of the front can be calculated 
similar to Eq. 1171 

Wi =^( =Us ). (29) 

71; 

The jump conditions can then be rewritten to give the fol- 
lowing relations: 

(30) 



U s . C -2 

n c = no— = n ( = — ) 



"0 

al 
'a? 



ctBnfao 



■-' n C 2 . 



= (- 



a 

ui — it s • 

U B 



(31) 
(32) 
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Classic 


Exact 


Grid 


SPHI 


iVINE 


n c (cm -3 ) 


147 


137 


169 


155 


138 ± 6 


rij (cm -3 ) 


0.734 


0.747 


0.748 


0.75 


0.743 ± 0.01 


Uj (km s — 1 ) 


3.48 


3.37 


3.36 


3.43 


3.34 ± 0.18 


Ml (km s — 1 ) 


3.24 


3.12 






3.13 ± 0.04 



100 r 



Table 1. Comparison of analytical and numerical results for the 
test including hydrodynamics and ionization. The iVINE data 
is obtained from the 2 million p article run, the erro r s give n are 
lcr. The grid d ata is taken from iLefloch &: Laz arcff (1994[), the 
SPHI data from lKessel-Devnet fc Burkertl fcOOOf) . The analytical 
values differ from the previous work due to a higher accuracy in 
our calculations. 



To compare dir ectly to previous results , we used the 
initial conditions bv lLefloch fc Lazarefj (| 19941 ). The density 
is no = 100cm -3 , the temperature is T co id = 100K. The flux 
increases linearly with time at a constant rate of dF/dt = 
5.07 x 10 -8 cm -2 s -2 , starting from zero. As before the re- 
combination parameter is set to qb = 2.7 x 10 _13 cm 3 s _1 
and the ionized temperature is Thot = 10 4 K. Refinement is 
included. 

The simulations are performed using the individual par- 
ticle time-steps of VINE. For the d etermination of t h e time - 
step we use the criteria given in IWetzstein et all (|2008t) . 
Here we will only review briefly the parameters used. We 
use a combined time-step criterion based on the change in 
acceleration and velocity of the particle with an accuracy 
parameter of r acc = 1.0. In addition, a CFL criterion is 
used with a tolerance parameter of tcfl = 0.3 and the 
modifications discussed in Section 12.31 We also use an ad- 
ditional time-step criterion based on the maximum allowed 
change of the smoothing length (see IWetzstein et all 120081 
for details) with an accuracy parameter of 71, = 0.15. VINE 
employs a variable and time-dependent smoothing length, 
the number of neighbours of each particle is on average 
Wneigh = 50, but variations of ±20 are allowed. The artificial 
viscosity of the SPH method is included in the standard way 
jCingold fc Monaghanll 19831 ) with the parameters a = 1 and 
P — 2 as implemented in VINE. 

We performed simulations with 1 and 2 million particles 
in a cubic simulation domain. The particles are distributed 
randomly and then left to relax according to Section [3.1.11 
The higher resolution compared to the test in Section 13.11 
is necessary to follow the hydrodynamical shock precisely. 
As in the tests before, assuming a fully ionized gas with 
r\ — 1.0 and thereby using a value of Thot, a = 10 4 K for the 
calculation of the sound-speed of the hot gas, does not corre- 
spond to the simulations (see equation I22|l . Instead a much 
better agreement can be achieved when the real tempera- 
ture of the gas in the simulations, Th ot , r = 9200K, is used 
(since r\ — 0.92 on average in the ionized region). For this 
more realistic temperature the simulations are in very good 
agreement and converge towards the analytic solution with 
increasing resolution (see Fig. [SJ . This can also be seen in 
Table Q] 



3.3 Performance 

To test the performance of the parallel iVINE code with 
increasing number of processors we choose the simulation 




Figure 5. Number density versus position for the steady prop- 
agation of an ionizing front. The dashed line shows the classic 
solution, obtained by using a value of T^ ot !l = 10 4 K for the hot 
gas. The solid line corresponds to the analytic solution for a more 
realistic value of T^ ot r = 9200K for the hot gas. Blue and green 
lines show the simulations at different resolutions. 



described in Section|4]at a later stage and compute one time- 
step on different numbers of CPUs. The parallel scaling of 
the various parts of the underl ying VINE code is discussed 
in detail in iNelson et al.l (|2008l ). For our test, we use a SGI 
Altix 3700 Bx2 supercomputer. In total, the ionization uses 
only a few percent of the total computational time. The 
precise values range from 2.32% on 2 CPUs to 2.70% on 16 
CPUs and 2.86% on 32 CPUs. 

When refinement is used, these values change to 8.52% 
on 2 CPUs and 8.73% on 32 CPUs. Although the ionization 
takes up relatively more time in this case, the difference 
in the calculation time between the number of CPUs gets 
smaller. This is to be expected, as the refinement is calcu- 
lated inside the bins and this part of the implementation is 
parallelized very efficiently (each bin is independent of the 
other bins). 

This test shows that the additional cost of our im- 
plementation of ionizing radiation in SPH is always much 
smaller than the cost for other implemented physics, like 
gravity and hydrodynamics. In particular, our new ray- 
tracing-scheme sh ows a substantial speedup comp ared to 

where 



up compa 



the algorithm by iKessel-Devnet fc Burkerfj 
the path-finding alone took up about 50% of the total com- 
putational time (lKessel-Devnet|[l999l ) . Another drawback of 
their approach is that for every particle the optical depth 
is calculated along a path towards the source until a parti- 
cle closer to the source with an already calculated optical 
depth is f ound. This is a highly seria l appr oach and thus the 
scheme o f lKessel-Devnet fc Burkertl (2000) does not lend it- 
self easily to an efficient parallelization. 



4 RADIATION DRIVEN IMPLOSION 

As a first application of iVINE we model the radiation driven 
implosion of an otherwise stable molecular cloud cor e . This 
approach is very similar to lKessel-Devnet fc Burkertl (|2003l ) 
but at ten times higher mass resolution. A marginally stable 
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Bonnor-Ebert sphere (BES) jBonnorl |l956) with a radial 
pressure profile defined by 



1 d 

r 2 ^ r 



2 dp 

dr 



-4nGp 



(33) 



is exposed to UV-radiation from a nearby source. The tem- 
perature of the sphere is T — 10K, the peak density is 
n max = 10 3 cm -3 , and the gas is initially at rest (i.e. no tur- 
bulence). The total mass contained in the sphere is 96M0 
and the radius is 1.6pc. We embed the sphere into cold gas 
(10K) with a constant density corresponding to the cutoff- 
density at the edge of the sphere. These simulations where 
performed with 2.2 x fO 6 particles resulting in a particle 
mass of 7.2 x 1O _5 M0. Self gravity is included. The cool- 
ing timescale (t coo \ < 0.3kyr) is much shorter than any 
other timescale involved in our simulations (e.g. the cross- 
ing timescale of the hot gas is that ~ 70kyr) . Thus, we treat 
the non-ionized gas with an isothermal equation of state 
(7 = f). The ionized gas is assigned a temperature accord- 
ing to Eq. Owith T hot = f0 4 K and T coM = fOK and then 
treated isothermally as well. 

The artificial viscosity and the criteria for the individual 
time-steps are the same as in Section [3.2.11 (a = 1, j3 = 2, 
Tacc = 1., tcfl = 0.3 and Th = 0.15). In addition, we use a 
multipole acceptance criterion (MAC) for the tree based cal- 
culati on of gravitational forces according to ISpringel et all 
l|200ll ) as implemented bv lWetzstein et all l|2008l ) with a tree 
accuracy parameter of 8 — 5 x 10 -4 . The implementation 
of the SPH smoothing kernel and the gravitational soften- 
ing length in VINE are equal at all times. The number of 
neighbours is set to n nc i g h = 50 ± 20. The hydrodynamical 
boundaries are periodic in the y- and z- direction, and open 
in the x-direction. This resembles the situation around a 
massive O-star where the material is allowed to move freely 
in the radial direction while at the sides similar material 
is existing. Gravitational forces are calculated by just tak- 
ing into account the self-gravity of the gas and no external 
or boundary effects. This is reasonable as the total simu- 
lation time (< 600kyr) is much shorter than the free-fall 
time (iff m 1.5Myr). In the simulations the radiation is im- 
pinging from the left hand side. We perform three different 
simulations, differentiated by the penetration length in the 
surrounding medium relative to the box size C = x s /xbox'- 



Simulation HF (high flux): 
F = 9.0 x 10 9 photons cm~ 2 s" 1 => C i 
Simulation IF (intermediate flux): 
F = 4.5 x 10 9 photons cm" 2 s _1 => C ~ 0.5 
Simulation LF (low flux): 
Jo = 9.0 x 10 8 photons cm" 



1.0 



= > C « 0.1 



This corresponds to the molecular cloud being placed inside 
(HF), at the border (IF) and outside (LF) of the Stromgren 
sphere. The evolution of the BES for all three cases is shown 
in Fig. [5] 

4.1 Dynamical evolution 

The general evolution of a simulation of this kind is as 
follows: As soon as the simulation starts, a R-type ion- 
ization front is driven into the medium. As it can be ex- 
pected from Section [3.1.11 the front reaches the Stromgren 
radius x s of the diffuse gas within a few recombination 



timescales (5t ICC ~ 5kyr). After a sound crossing timescale 
(thot ~ 70kyr) the hot gas reacts to its change in pressure 
and starts to drive a shock front into the cold gas - a D-type 
front evolves. This front will affect the morphology of the 
BES. In the following we describe the individual cases in 
more detail. 



4-1.1 Simulation HF (high flux) 

Due to the high flux (see Fig. [6] first column), the R-type 
front is able to propagate very far into the simulation do- 
main. A bow-like shock structure around the edge of the 
BES evolves (t « lOOkyr) . The shock front running into the 
denser parts of the cloud is slowed down, so that the front 
starts to "wrap around" the cloud. Soon the two flanks are 
approaching each other while the center of the shock is held 
back by the dense innermost region (Fig. [6] third row, first 
column, t ~ lOOkyr). As the two sides finally collide an 
elongated filament forms which is gravitationally unstable. 
In Fig. [7] we show this filament at the final stage of our sim- 
ulations. In comparison runs without self-gravity the two 
shock fronts cross each other and the cloud disperses. With 
the inclusion of self-gravity the filament becomes gravita- 
tional unstable and is triggered into collapse. In fact the 
core fragments into several objects, as will be discussed 
in a subsequent paper . The resolution limit according to 
iBate fc Burkertl (|l997h is nmax — 2 x 10 10 cm 3 in these 
simulations. As soon as this limit is reached the local Jeans 
mass is smaller than the mass of 100 particles and artifical 
fragmentation can occur. Thus, we stop the simulations at 
this point. 



4-1.2 Simulation IF (intermediate flux) 

With an intermediate flux the R-type front penetrates much 
less into the gaseous medium (see Fig. [5] second column). 
Thus, the front does not wrap around the sphere as much 
as in Simulation HF. As soon as the hot gas reacts to its 
increased temperature a flattened shock with a much smaller 
curvature than in Simulation HF forms. In addition, the 
motion in the hot gas forces the flanks on both sides of the 
main shock inwards, which can be seen in the velocities of 
the hot gas in Fig. [6] (third row, second column, t m lOOkyr). 
These motions are due to the periodic boundaries on the 
upper and lower edge. Otherwise the gas could stream away 
freely. However these boundaries are justified by the fact 
that the molecular cloud is completely surrounding the O- 
star. In the further evolution the flanks approach each other 
similar to Simulation HF and the central region becomes 
unstable and fragments (see Fig. 0). 



4-1-3 Simulation LF (low flux) 

The very low flux in this case only leads to a R-front which 
barely reaches the sphere (see Fig. [6] third column). There- 
fore, the D-front starts as a nearly plane-parallel shock wave 
in front of the BES. This shock sweeps away much more ma- 
terial than in the high and intermediate flux cases, where 
the material is concentrated in the center. Nevertheless, as 
the shock propagates further the very center of the sphere 
gets compressed and becomes gravitationally unstable. In 



iVINE - Ionization in the tree/SPH code VINE 9 



High flux (HF) Intermediate flux (IF) Low flux (LF) 




-2.0 -1.0 0.0 1.0 2.0 -2.0 -1.0 0.0 1.0 2.0 -2.0 -1.0 0.0 1.0 2.0 

x [pc] x [pc] x [pc] 



Figure 6. Time evolution of the driven collapse of a Bonnor-Ebert sphere ionized by a source with high flux (first column), intermediate 
flux (second column) and low flux (third column). The simulation volume is a cube with sides 4pc long, the UV-radiation is impinging 
from the left hand side. Color coded is the density of the central cold gas slab. Yellow arrows denote the velocities of the hot gas, black 
arrows the motion of the cold gas. Density and velocities are averaged across a slice of 0.125pc in the z-direction. Each row shows the 



10 M. Gritschneder, T. Naab, A. Burkert, S. Walch, F. Heitsch, M. Wetzstein 



contrast to Simulations HF and IF there is no sign of frag- 
mentation in the unstable region. 

4.2 Structure, collapse timescales and final mass 
assembled 

A close look at the final structure of the collapsing fil- 
aments of the three simulations (Fig. [TJl reveals that in 
all three simulations the core forms at the tip of an elon- 
gated filament, which might be eroded in the future. This 
matches exactly the observed head to tail structure de- 
scribed in Section [1] The core regions have an extent of 
just 0.02 — 0.05pc, which cor r espon ds very closely to e.g. the 
findings of iMotte fc Andrei ((2001) in the the Perseus star 
cluster. They observe 8 Class protostars with compact 
envelopes (R ou t < 10 4 AU w 0.05pc). In addition they are 
denser by a factor of 3-12 than it would be expected from 
the standard collapse model, which would suggest densities 
of n sa 10 6 cm~ 3 (see e .g. Walch et al. 2008, in preparation). 
IMotte fc Andrei l|200ll ) suggest that this higher densities are 
due to external disturbances initiating the collapse, which 
agrees very well with our simulations. Following the obser- 
vations we define a core as all material with a density higher 
than n C rit = 10 7 cm~ 3 in a region of i? cr it = 0.02pc (which 
is roughly a Jeans length at a density of n C rit) around the 
peak density. 

We plot the evolution of the maximum density in Fig. 
[S] In all three simulations after the first phase of compres- 
sion by the hot gas a meta-stable phase at densities between 
10 6 cm -3 and 10 7 cm -3 can be seen. This fits very nicely to 
the structure of observed cores described above. The dura- 
tion of this phase depends on the initial flux (HF: 90kyr, 
IF: 155kyr, LF: 290kyr). In addition, we find evidence that 
the filaments collapse earlier in cases with a higher flux. The 
collapse happens at t — 200kyr, t = 280kyr and t = 600kyr 
in Simulation HF, IF and LF, respectively. Observations of 
triggered star f ormation tend to show the same trend (see 
iLee et al.|[2005l . Ilkeda et al.ll2008l ) - the younger the star, the 
further it is away from the ionizing source. This can not be 
explained by just attributing it to the speed of the R-type 
front. As seen in Section [3.1. li the crossing time for the R- 
type front is of the order of a few kyr, whereas any observed 
age-spread is of the order of several hundred kyr. To explain 
this huge spread the position of the density enhancement 
relative to the Stromgren radius has to be considered. As 
we show decreasing the flux and thereby increasing the dis- 
tance to the source can delay collapse and star formation by 
0.08 - 0.4Myr. 

In IC 1396N iGetman et ail (|2007l ) report a T Tauri 
(Class II and Class III stars) population with ages ~ 0.5 — 1 
Myr. In addition, 0.3 — 0.5pc further away from the ionizing 
source HD 206267 (an 06.5f-type star), there is an embed- 
ded population of Class 0/1 protostars with ages « O.lMyr. 
This can be compared to our simulations where e.g Simu- 
lation IF represents gas clumps closer to the source which 
start to form stars 0.3Myr earlier than Simulation LF. So 
at the time the embedded stars start to form in Simulation 
LF the stars of Simulation IF would be no longer embedded 
and represent the Class II/III stars population. In fact in 
our simulations the spread of a few hundred kyr is smaller 
than in the observations. This difference could be attributed 
to the classification of the protostars as discussed below. 



High flux (HF) 




Intermediate flux (IF) 




Low flux (LF) 




0.5 1.0 1.5 a.o 2.5 

x [pc] 



Figure 7. Final stage of the three simulations. Color coded is the 
density in the central slab. Yellow arrows denote the velocities of 
the hot gas, black arrows the motion of the cold gas. Density and 
velocities are averaged across a slice of 0.0625pc in the z-direction. 
The time of the collapse as well as the displacement of the frag- 
ment clearly depend on the initial flux. Furthermore the velocity 
of the cold gas (black arrows) is decreasing with decreasing flux. 
The core always forms at the very tip of the filament. 



Besides the age-spread one can look at the velocities 
of the front and core. From Fig. [Tj] it can be seen that the 
shock front travels with a speed of 3 — 7km/s, depending 
on the initial flux. M ost observational estim ates provide a 
front speed < lkm/s l|Thompson et al1l2004l ). leading to a 
difference of almost an order of magnitude between observa- 
tions and simulation s which has been noted before (see e.g. 
IGetman et ail 120071 ). The age estimates of the YSOs are 
mainly based on their classification by the spectral energy 
distribution (SED). A Class 0/1 object is deeply embedded, 
therefore the short micrometer wavelengths are much weaker 
due to absorptions in the envelope when compared to Class 
II/III objects. This allows for a clear distinction between 
both types e.g. in the I RAC [3.6] -[4.5] versus [5.8] - [8] color 
diagram as shown by lHartmann et al.l (|2005l ) . We suggest 
that in the case of triggered star formation, the ionizing ra- 
diation could strip the envelope of a YSO, thereby unveiling 
the central object in short micrometer wavelengths. Thus, 
the observed Class II/III SED could be caused by a much 
younger Class 0/1 protostar with a removed envelope. This 
would reduce the estimated age spread, thereby increasing 
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Figure 8. The maximum number density versus time for the 
three different simulations. In the higher flux cases HF and IF 
the collapse happens much earlier than in the low flux case LF. 
The dash-dotted line represents the resolution limit as given by 
iBate fc Burked lll997f) . 

the estimated speed of the shock front and finally leading 
observations and simulations to agree. This assumption will 
be subject to further examination. 

A dependence on distance is also seen for the veloci- 
ties of the cold filaments (see the velocities of the cold gas 
(black arrows) in Fig. [7J . The precise velocities of the cores 
in the Simulations HF, IF and LF are 8.4km/s, 7.6km/s 
and 5.1km/s, respectively. Again, the closer the core is to 
the OB-association the higher is its velocity. Although this 
small differences are not observable yet it is worth notic- 
ing that the cores themselves have bulk velocities which are 
slightly higher (by about 0.5 — lkm/s) than the rest of the 
filament. However, this effect may get weaker as the core 
gets slowed down while sweeping up the rest of the filament. 

The final mass assembled does not show a dependence 
on the initial distance. In Simulation HF the core consists 
of 6.0M Q in Simulation IF of 7.4M and in Simulation LF 
of 2.8Mq. The filaments as a whole have masses of 61.5Mq, 
75.3Mq or 67.4M0, respectively. It is obvious that the most 
effective scenario is Simulation IF. Here, the ionization en- 
compasses most of the sphere and thus the shock front is 
not nearly as plane-parallel as in Simulation LF and does 
not sweep away as much material. On the other hand, less 
material gets evaporated by the ionization since the flux is 
lower than in Simulation HF. Overall the final masses of the 
collapsing cores fit the observations well. Assuming a star 
formation efficiency of 30% (see e.g. lLada et al.ll200cj ). we 
find masses from 0.84Mq to 2.2Mq which agrees with the 
observed ran ge from classical T Tauri up to Herb ig Ae/Be 
stars (see e.s JLee fc Chenll2007l. [Snider et al.ll2007h . 



5 SUMMARY AND DISCUSSION 

We present iVINE, a new implementation of UV-radiation 
into the tree-SPH code VINE. It uses a plane-parallel geom- 
etry which renders the code most suitable to perform high 
resolution studies of the small scale effects of e.g. ionization 
and turbulence in the surrounding of young massive stars. 



It is efficiently parallelized and very fast, as only 2%-8% of 
the total computational time are used for the calculation 
of ionization. The comparison with analytic solutions shows 
that iVINE treats time-dependent ionization as well as the 
resulting heating effects precisely and convergently. 

We base our numerical implementation of ionizing ra- 
diation on several assumptions. First, we use a simplified 
prescription for the radiative transfer by e.g assuming a 
monochromatic flux. Second, we neglect UV absorption by 
dust, which would lower the total UV flux. Third, we do not 
include a full treatment of recombination zones. In our sim- 
ulations the ionized gas which gets shaded is assumed to re- 
combine immediately. In addition, the gas in the shaded re- 
gions does not get heated by irradiation from the hot gas sur- 
rounding it. These effects require a precise time-dependent 
treatment of heating and cooling processes by ionization and 
recombination as well as a treatment for the scattering of 
photons. An implementation of this effects is planned in a 
future version of the code. 

As an application we investigate radiation driven implo- 
sion of a marginally stable Bonnor-Ebert sphere. We show 
that these spheres are indeed driven into gravitational col- 
lapse. The resulting cores are in the observed mass range. 
They are more compact and a factor of ~ 10 more dense 
than it would be expected in a more quiescent environment. 
This fact fits very well with the observations of star forma- 
tion in a clustered environment. By comparing simulations 
with three different UV-fluxes we show that there is a clear 
dependence of the final mass and the age of the collapsed 
core on the position of the preexisting density enhancement 
relative to the Stromgren radius. Our findings that the on- 
set of star formation is delayed by 0.08 — 0.4Myr, depending 
on the position, are in good agreement with observations of 
the age spread in bright rimmed clouds. The velocity of the 
triggering shock is an order of magnitude higher than the 
observational estimates. This discrepancy has been noted 
before. We suggest that this can be attributed to the ioniz- 
ing radiation stripping the envelope from a Class 0/1 star. 
Thereby it might be classified as an Class II/III star, leading 
to an higher age-spread between the observed objects. Cor- 
recting for this effect would increase the estimated velocity 
of the shock front and thus lead simulations and observa- 
tions towards agreement. 
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